!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief 2- and 3-center electron repulsion integral routines based on libint2
!>        Currently available operators: Coulomb, Truncated Coulomb, Short Range (erfc), Overlap
!> \author A. Bussy (05.2019)
! **************************************************************************************************

MODULE libint_2c_3c
   USE gamma,                           ONLY: fgamma => fgamma_0
   USE input_constants,                 ONLY: do_potential_coulomb,&
                                              do_potential_id,&
                                              do_potential_long,&
                                              do_potential_mix_cl_trunc,&
                                              do_potential_short,&
                                              do_potential_truncated
   USE kinds,                           ONLY: default_path_length,&
                                              dp
   USE libint_wrapper,                  ONLY: cp_libint_get_2eri_derivs,&
                                              cp_libint_get_2eris,&
                                              cp_libint_get_3eri_derivs,&
                                              cp_libint_get_3eris,&
                                              cp_libint_set_params_eri,&
                                              cp_libint_set_params_eri_deriv,&
                                              cp_libint_t,&
                                              prim_data_f_size
   USE mathconstants,                   ONLY: pi
   USE orbital_pointers,                ONLY: nco,&
                                              ncoset
   USE t_c_g0,                          ONLY: get_lmax_init,&
                                              t_c_g0_n
#include "./base/base_uses.f90"

   IMPLICIT NONE
   PRIVATE

   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libint_2c_3c'

   PUBLIC :: eri_3center, eri_2center, cutoff_screen_factor, libint_potential_type, &
             eri_3center_derivs, eri_2center_derivs, compare_potential_types

   ! For screening of integrals with a truncated potential, it is important to use a slightly larger
   ! cutoff radius due to the discontinuity of the truncated Coulomb potential at the cutoff radius.
   REAL(KIND=dp), PARAMETER :: cutoff_screen_factor = 1.0001_dp

   TYPE :: params_2c
      INTEGER                               :: m_max = 0
      REAL(dp)                              :: ZetaInv = 0.0_dp, EtaInv = 0.0_dp, ZetapEtaInv = 0.0_dp, Rho = 0.0_dp
      REAL(dp), DIMENSION(3)                :: W = 0.0_dp
      REAL(dp), DIMENSION(prim_data_f_size) :: Fm = 0.0_dp
   END TYPE params_2c

   TYPE :: params_3c
      INTEGER                               :: m_max = 0
      REAL(dp)                              :: ZetaInv = 0.0_dp, EtaInv = 0.0_dp, ZetapEtaInv = 0.0_dp, Rho = 0.0_dp
      REAL(dp), DIMENSION(3)                :: Q = 0.0_dp, W = 0.0_dp
      REAL(dp), DIMENSION(prim_data_f_size) :: Fm = 0.0_dp
   END TYPE params_3c

   ! A potential type based on the hfx_potential_type concept, but only for implemented
   ! 2- and 3-center operators
   TYPE :: libint_potential_type
      INTEGER                               :: potential_type = do_potential_coulomb
      REAL(dp)                              :: omega = 0.0_dp !SR: erfc(w*r)/r
      REAL(dp)                              :: cutoff_radius = 0.0_dp !TC cutoff/effective SR range
      CHARACTER(default_path_length)        :: filename = ""
      REAL(dp)                              :: scale_coulomb = 0.0_dp ! Only, for WFC methods
      REAL(dp)                              :: scale_longrange = 0.0_dp ! Only, for WFC methods
   END TYPE libint_potential_type

CONTAINS

! **************************************************************************************************
!> \brief Computes the 3-center electron repulsion integrals (ab|c) for a given set of cartesian
!>        gaussian orbitals
!> \param int_abc the integrals as array of cartesian orbitals (allocated before hand)
!> \param la_min ...
!> \param la_max ...
!> \param npgfa ...
!> \param zeta ...
!> \param rpgfa ...
!> \param ra ...
!> \param lb_min ...
!> \param lb_max ...
!> \param npgfb ...
!> \param zetb ...
!> \param rpgfb ...
!> \param rb ...
!> \param lc_min ...
!> \param lc_max ...
!> \param npgfc ...
!> \param zetc ...
!> \param rpgfc ...
!> \param rc ...
!> \param dab ...
!> \param dac ...
!> \param dbc ...
!> \param lib the libint_t object for evaluation (assume that it is initialized outside)
!> \param potential_parameter the info about the potential
!> \param int_abc_ext the extremal value of int_abc, i.e., MAXVAL(ABS(int_abc))
!> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
!>       the libint library must be static initialized, and in case of truncated Coulomb operator,
!>       the latter must be initialized too
! **************************************************************************************************
   SUBROUTINE eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, &
                          lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
                          lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
                          dab, dac, dbc, lib, potential_parameter, &
                          int_abc_ext)

      REAL(dp), DIMENSION(:, :, :), INTENT(INOUT)        :: int_abc
      INTEGER, INTENT(IN)                                :: la_min, la_max, npgfa
      REAL(dp), DIMENSION(:), INTENT(IN)                 :: zeta, rpgfa
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: ra
      INTEGER, INTENT(IN)                                :: lb_min, lb_max, npgfb
      REAL(dp), DIMENSION(:), INTENT(IN)                 :: zetb, rpgfb
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: rb
      INTEGER, INTENT(IN)                                :: lc_min, lc_max, npgfc
      REAL(dp), DIMENSION(:), INTENT(IN)                 :: zetc, rpgfc
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: rc
      REAL(KIND=dp), INTENT(IN)                          :: dab, dac, dbc
      TYPE(cp_libint_t), INTENT(INOUT)                   :: lib
      TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter
      REAL(dp), INTENT(INOUT), OPTIONAL                  :: int_abc_ext

      INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, ipgf, j, &
         jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
      REAL(dp)                                           :: dr_ab, dr_ac, dr_bc, zeti, zetj, zetk
      REAL(dp), DIMENSION(:), POINTER                    :: p_work
      TYPE(params_3c), POINTER                           :: params

      NULLIFY (params, p_work)
      ALLOCATE (params)

      dr_ab = 0.0_dp
      dr_bc = 0.0_dp
      dr_ac = 0.0_dp

      op = potential_parameter%potential_type

      IF (op == do_potential_truncated .OR. op == do_potential_short &
          .OR. op == do_potential_mix_cl_trunc) THEN
         dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
         dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
      ELSEIF (op == do_potential_coulomb) THEN
         dr_bc = 1000000.0_dp
         dr_ac = 1000000.0_dp
      END IF

      IF (PRESENT(int_abc_ext)) THEN
         int_abc_ext = 0.0_dp
      END IF

      !Note: we want to compute all possible integrals based on the 3-centers (ab|c) before
      !      having to switch to (ba|c) (or the other way around) due to angular momenta in libint
      !      For a triplet of centers (k|ji), we can only compute integrals for which lj >= li

      !Looping over the pgfs
      DO ipgf = 1, npgfa
         zeti = zeta(ipgf)
         a_start = (ipgf - 1)*ncoset(la_max)

         DO jpgf = 1, npgfb

            ! screening
            IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE

            zetj = zetb(jpgf)
            b_start = (jpgf - 1)*ncoset(lb_max)

            DO kpgf = 1, npgfc

               ! screening
               IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
               IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE

               zetk = zetc(kpgf)
               c_start = (kpgf - 1)*ncoset(lc_max)

               !start with all the (c|ba) integrals (standard order) and keep to lb >= la
               CALL set_params_3c(lib, ra, rb, rc, zeti, zetj, zetk, la_max, lb_max, lc_max, &
                                  potential_parameter=potential_parameter, params_out=params)

               DO li = la_min, la_max
                  a_offset = a_start + ncoset(li - 1)
                  ncoa = nco(li)
                  DO lj = MAX(li, lb_min), lb_max
                     b_offset = b_start + ncoset(lj - 1)
                     ncob = nco(lj)
                     DO lk = lc_min, lc_max
                        c_offset = c_start + ncoset(lk - 1)
                        ncoc = nco(lk)

                        a_mysize(1) = ncoa*ncob*ncoc
                        CALL cp_libint_get_3eris(li, lj, lk, lib, p_work, a_mysize)

                        IF (PRESENT(int_abc_ext)) THEN
                           DO k = 1, ncoc
                              p1 = (k - 1)*ncob
                              DO j = 1, ncob
                                 p2 = (p1 + j - 1)*ncoa
                                 DO i = 1, ncoa
                                    p3 = p2 + i
                                    int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
                                    int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
                                 END DO
                              END DO
                           END DO
                        ELSE
                           DO k = 1, ncoc
                              p1 = (k - 1)*ncob
                              DO j = 1, ncob
                                 p2 = (p1 + j - 1)*ncoa
                                 DO i = 1, ncoa
                                    p3 = p2 + i
                                    int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
                                 END DO
                              END DO
                           END DO
                        END IF

                     END DO !lk
                  END DO !lj
               END DO !li

               !swap centers 3 and 4 to compute (c|ab) with lb < la
               CALL set_params_3c(lib, rb, ra, rc, params_in=params)

               DO lj = lb_min, lb_max
                  b_offset = b_start + ncoset(lj - 1)
                  ncob = nco(lj)
                  DO li = MAX(lj + 1, la_min), la_max
                     a_offset = a_start + ncoset(li - 1)
                     ncoa = nco(li)
                     DO lk = lc_min, lc_max
                        c_offset = c_start + ncoset(lk - 1)
                        ncoc = nco(lk)

                        a_mysize(1) = ncoa*ncob*ncoc
                        CALL cp_libint_get_3eris(lj, li, lk, lib, p_work, a_mysize)

                        IF (PRESENT(int_abc_ext)) THEN
                           DO k = 1, ncoc
                              p1 = (k - 1)*ncoa
                              DO i = 1, ncoa
                                 p2 = (p1 + i - 1)*ncob
                                 DO j = 1, ncob
                                    p3 = p2 + j
                                    int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
                                    int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
                                 END DO
                              END DO
                           END DO
                        ELSE
                           DO k = 1, ncoc
                              p1 = (k - 1)*ncoa
                              DO i = 1, ncoa
                                 p2 = (p1 + i - 1)*ncob
                                 DO j = 1, ncob
                                    p3 = p2 + j
                                    int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
                                 END DO
                              END DO
                           END DO
                        END IF

                     END DO !lk
                  END DO !li
               END DO !lj

            END DO !kpgf
         END DO !jpgf
      END DO !ipgf

      DEALLOCATE (params)

   END SUBROUTINE eri_3center

! **************************************************************************************************
!> \brief Sets the internals of the cp_libint_t object for integrals of type (k|ji)
!> \param lib ..
!> \param ri ...
!> \param rj ...
!> \param rk ...
!> \param zeti ...
!> \param zetj ...
!> \param zetk ...
!> \param li_max ...
!> \param lj_max ...
!> \param lk_max ...
!> \param potential_parameter ...
!> \param params_in external parameters to use for libint
!> \param params_out returns the libint parameters computed based on the other arguments
!> \note The use of params_in and params_out comes from the fact that one might have to swap
!>       centers 3 and 4 because of angular momenta and pretty much all the parameters of libint
!>       remain the same upon such a change => might avoid recomputing things over and over again
! **************************************************************************************************
   SUBROUTINE set_params_3c(lib, ri, rj, rk, zeti, zetj, zetk, li_max, lj_max, lk_max, &
                            potential_parameter, params_in, params_out)

      TYPE(cp_libint_t), INTENT(INOUT)                   :: lib
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: ri, rj, rk
      REAL(dp), INTENT(IN), OPTIONAL                     :: zeti, zetj, zetk
      INTEGER, INTENT(IN), OPTIONAL                      :: li_max, lj_max, lk_max
      TYPE(libint_potential_type), INTENT(IN), OPTIONAL  :: potential_parameter
      TYPE(params_3c), OPTIONAL, POINTER                 :: params_in, params_out

      INTEGER                                            :: l
      LOGICAL                                            :: use_gamma
      REAL(dp)                                           :: gammaq, omega2, omega_corr, omega_corr2, &
                                                            prefac, R, S1234, T, tmp
      REAL(dp), ALLOCATABLE, DIMENSION(:)                :: Fm
      TYPE(params_3c), POINTER                           :: params

      !Assume that one of params_in or params_out is present, and that in the latter case, all
      !other optional arguments are here

      !The internal structure of libint2 is based on 4-center integrals
      !For 3-center, one of those is a dummy center
      !The integral is assumed to be (k|ji) where the centers are ordered as:
      !k -> 1, j -> 3 and i -> 4 (the center #2 is the dummy center)

      !If external parameters are given, just use them
      IF (PRESENT(params_in)) THEN
         params => params_in

         !If no external parameters to use, compute them
      ELSE
         params => params_out

         !Note: some variable of 4-center integrals simplify with a dummy center:
         !      P -> rk, gammap -> zetk
         params%m_max = li_max + lj_max + lk_max
         gammaq = zeti + zetj
         params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
         params%ZetapEtaInv = 1._dp/(zetk + gammaq)

         params%Q = (zeti*ri + zetj*rj)*params%EtaInv
         params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
         params%Rho = zetk*gammaq/(zetk + gammaq)

         params%Fm = 0.0_dp
         SELECT CASE (potential_parameter%potential_type)
         CASE (do_potential_coulomb)
            T = params%Rho*SUM((params%Q - rk)**2)
            S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
            prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234

            CALL fgamma(params%m_max, T, params%Fm)
            params%Fm = prefac*params%Fm
         CASE (do_potential_truncated)
            R = potential_parameter%cutoff_radius*SQRT(params%Rho)
            T = params%Rho*SUM((params%Q - rk)**2)
            S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
            prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234

            CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
            CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
            IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
            params%Fm = prefac*params%Fm
         CASE (do_potential_short)
            T = params%Rho*SUM((params%Q - rk)**2)
            S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
            prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234

            CALL fgamma(params%m_max, T, params%Fm)

            omega2 = potential_parameter%omega**2
            omega_corr2 = omega2/(omega2 + params%Rho)
            omega_corr = SQRT(omega_corr2)
            T = T*omega_corr2
            ALLOCATE (Fm(prim_data_f_size))

            CALL fgamma(params%m_max, T, Fm)
            tmp = -omega_corr
            DO l = 1, params%m_max + 1
               params%Fm(l) = params%Fm(l) + Fm(l)*tmp
               tmp = tmp*omega_corr2
            END DO
            params%Fm = prefac*params%Fm
         CASE (do_potential_mix_cl_trunc)
            R = potential_parameter%cutoff_radius*SQRT(params%Rho)
            T = params%Rho*SUM((params%Q - rk)**2)
            S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
            prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234

            CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
            CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
            IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)

            ALLOCATE (Fm(prim_data_f_size))
            CALL fgamma(params%m_max, T, Fm)
            DO l = 1, params%m_max + 1
               params%Fm(l) = params%Fm(l) &
                              *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
                              - Fm(l)*potential_parameter%scale_longrange
            END DO
            DEALLOCATE (Fm)

            omega2 = potential_parameter%omega**2
            omega_corr2 = omega2/(omega2 + params%Rho)
            omega_corr = SQRT(omega_corr2)
            T = T*omega_corr2

            ALLOCATE (Fm(prim_data_f_size))
            CALL fgamma(params%m_max, T, Fm)
            tmp = omega_corr
            DO l = 1, params%m_max + 1
               params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
               tmp = tmp*omega_corr2
            END DO
            params%Fm = prefac*params%Fm
         CASE (do_potential_id)
            S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
                        - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
            prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234

            params%Fm(:) = prefac
         CASE DEFAULT
            CPABORT("Requested operator NYI")
         END SELECT

      END IF

      CALL cp_libint_set_params_eri(lib, rk, rk, rj, ri, params%ZetaInv, params%EtaInv, &
                                    params%ZetapEtaInv, params%Rho, rk, params%Q, params%W, &
                                    params%m_max, params%Fm)

   END SUBROUTINE set_params_3c

! **************************************************************************************************
!> \brief Computes the derivatives of the 3-center electron repulsion integrals (ab|c) for a given
!>        set of cartesian gaussian orbitals. Returns x,y,z derivatives for 1st and 2nd center
!> \param der_abc_1 the derivatives for the 1st center (allocated before hand)
!> \param der_abc_2 the derivatives for the 2nd center (allocated before hand)
!> \param la_min ...
!> \param la_max ...
!> \param npgfa ...
!> \param zeta ...
!> \param rpgfa ...
!> \param ra ...
!> \param lb_min ...
!> \param lb_max ...
!> \param npgfb ...
!> \param zetb ...
!> \param rpgfb ...
!> \param rb ...
!> \param lc_min ...
!> \param lc_max ...
!> \param npgfc ...
!> \param zetc ...
!> \param rpgfc ...
!> \param rc ...
!> \param dab ...
!> \param dac ...
!> \param dbc ...
!> \param lib the libint_t object for evaluation (assume that it is initialized outside)
!> \param potential_parameter the info about the potential
!> \param der_abc_1_ext the extremal value of der_abc_1, i.e., MAXVAL(ABS(der_abc_1))
!> \param der_abc_2_ext ...
!> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
!>       the libint library must be static initialized, and in case of truncated Coulomb operator,
!>       the latter must be initialized too. Note that the derivative wrt to the third center
!>       can be obtained via translational invariance
! **************************************************************************************************
   SUBROUTINE eri_3center_derivs(der_abc_1, der_abc_2, &
                                 la_min, la_max, npgfa, zeta, rpgfa, ra, &
                                 lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
                                 lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
                                 dab, dac, dbc, lib, potential_parameter, &
                                 der_abc_1_ext, der_abc_2_ext)

      REAL(dp), DIMENSION(:, :, :, :), INTENT(INOUT)     :: der_abc_1, der_abc_2
      INTEGER, INTENT(IN)                                :: la_min, la_max, npgfa
      REAL(dp), DIMENSION(:), INTENT(IN)                 :: zeta, rpgfa
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: ra
      INTEGER, INTENT(IN)                                :: lb_min, lb_max, npgfb
      REAL(dp), DIMENSION(:), INTENT(IN)                 :: zetb, rpgfb
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: rb
      INTEGER, INTENT(IN)                                :: lc_min, lc_max, npgfc
      REAL(dp), DIMENSION(:), INTENT(IN)                 :: zetc, rpgfc
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: rc
      REAL(KIND=dp), INTENT(IN)                          :: dab, dac, dbc
      TYPE(cp_libint_t), INTENT(INOUT)                   :: lib
      TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter
      REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL      :: der_abc_1_ext, der_abc_2_ext

      INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, i_deriv, &
         ipgf, j, jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
      INTEGER, DIMENSION(3)                              :: permute_1, permute_2
      LOGICAL                                            :: do_ext
      REAL(dp)                                           :: dr_ab, dr_ac, dr_bc, zeti, zetj, zetk
      REAL(dp), DIMENSION(3)                             :: der_abc_1_ext_prv, der_abc_2_ext_prv
      REAL(dp), DIMENSION(:, :), POINTER                 :: p_deriv
      TYPE(params_3c), POINTER                           :: params

      NULLIFY (params, p_deriv)
      ALLOCATE (params)

      permute_1 = [4, 5, 6]
      permute_2 = [7, 8, 9]

      dr_ab = 0.0_dp
      dr_bc = 0.0_dp
      dr_ac = 0.0_dp

      op = potential_parameter%potential_type

      IF (op == do_potential_truncated .OR. op == do_potential_short &
          .OR. op == do_potential_mix_cl_trunc) THEN
         dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
         dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
      ELSEIF (op == do_potential_coulomb) THEN
         dr_bc = 1000000.0_dp
         dr_ac = 1000000.0_dp
      END IF

      do_ext = .FALSE.
      IF (PRESENT(der_abc_1_ext) .OR. PRESENT(der_abc_2_ext)) do_ext = .TRUE.
      der_abc_1_ext_prv = 0.0_dp
      der_abc_2_ext_prv = 0.0_dp

      !Note: we want to compute all possible integrals based on the 3-centers (ab|c) before
      !      having to switch to (ba|c) (or the other way around) due to angular momenta in libint
      !      For a triplet of centers (k|ji), we can only compute integrals for which lj >= li

      !Looping over the pgfs
      DO ipgf = 1, npgfa
         zeti = zeta(ipgf)
         a_start = (ipgf - 1)*ncoset(la_max)

         DO jpgf = 1, npgfb

            ! screening
            IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE

            zetj = zetb(jpgf)
            b_start = (jpgf - 1)*ncoset(lb_max)

            DO kpgf = 1, npgfc

               ! screening
               IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
               IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE

               zetk = zetc(kpgf)
               c_start = (kpgf - 1)*ncoset(lc_max)

               !start with all the (c|ba) integrals (standard order) and keep to lb >= la
               CALL set_params_3c_deriv(lib, ra, rb, rc, zeti, zetj, zetk, la_max, lb_max, lc_max, &
                                        potential_parameter=potential_parameter, params_out=params)

               DO li = la_min, la_max
                  a_offset = a_start + ncoset(li - 1)
                  ncoa = nco(li)
                  DO lj = MAX(li, lb_min), lb_max
                     b_offset = b_start + ncoset(lj - 1)
                     ncob = nco(lj)
                     DO lk = lc_min, lc_max
                        c_offset = c_start + ncoset(lk - 1)
                        ncoc = nco(lk)

                        a_mysize(1) = ncoa*ncob*ncoc

                        CALL cp_libint_get_3eri_derivs(li, lj, lk, lib, p_deriv, a_mysize)

                        IF (do_ext) THEN
                           DO i_deriv = 1, 3
                              DO k = 1, ncoc
                                 p1 = (k - 1)*ncob
                                 DO j = 1, ncob
                                    p2 = (p1 + j - 1)*ncoa
                                    DO i = 1, ncoa
                                       p3 = p2 + i

                                       der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
                                          p_deriv(p3, permute_2(i_deriv))
                                       der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
                                                                        ABS(p_deriv(p3, permute_2(i_deriv))))

                                       der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
                                          p_deriv(p3, permute_1(i_deriv))
                                       der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
                                                                        ABS(p_deriv(p3, permute_1(i_deriv))))

                                    END DO
                                 END DO
                              END DO
                           END DO
                        ELSE
                           DO i_deriv = 1, 3
                              DO k = 1, ncoc
                                 p1 = (k - 1)*ncob
                                 DO j = 1, ncob
                                    p2 = (p1 + j - 1)*ncoa
                                    DO i = 1, ncoa
                                       p3 = p2 + i

                                       der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
                                          p_deriv(p3, permute_2(i_deriv))

                                       der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
                                          p_deriv(p3, permute_1(i_deriv))
                                    END DO
                                 END DO
                              END DO
                           END DO
                        END IF

                        DEALLOCATE (p_deriv)
                     END DO !lk
                  END DO !lj
               END DO !li

               !swap centers 3 and 4 to compute (c|ab) with lb < la
               CALL set_params_3c_deriv(lib, rb, ra, rc, zetj, zeti, zetk, params_in=params)

               DO lj = lb_min, lb_max
                  b_offset = b_start + ncoset(lj - 1)
                  ncob = nco(lj)
                  DO li = MAX(lj + 1, la_min), la_max
                     a_offset = a_start + ncoset(li - 1)
                     ncoa = nco(li)
                     DO lk = lc_min, lc_max
                        c_offset = c_start + ncoset(lk - 1)
                        ncoc = nco(lk)

                        a_mysize(1) = ncoa*ncob*ncoc
                        CALL cp_libint_get_3eri_derivs(lj, li, lk, lib, p_deriv, a_mysize)

                        IF (do_ext) THEN
                           DO i_deriv = 1, 3
                              DO k = 1, ncoc
                                 p1 = (k - 1)*ncoa
                                 DO i = 1, ncoa
                                    p2 = (p1 + i - 1)*ncob
                                    DO j = 1, ncob
                                       p3 = p2 + j

                                       der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
                                          p_deriv(p3, permute_1(i_deriv))

                                       der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
                                                                        ABS(p_deriv(p3, permute_1(i_deriv))))

                                       der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
                                          p_deriv(p3, permute_2(i_deriv))

                                       der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
                                                                        ABS(p_deriv(p3, permute_2(i_deriv))))
                                    END DO
                                 END DO
                              END DO
                           END DO
                        ELSE
                           DO i_deriv = 1, 3
                              DO k = 1, ncoc
                                 p1 = (k - 1)*ncoa
                                 DO i = 1, ncoa
                                    p2 = (p1 + i - 1)*ncob
                                    DO j = 1, ncob
                                       p3 = p2 + j

                                       der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
                                          p_deriv(p3, permute_1(i_deriv))

                                       der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
                                          p_deriv(p3, permute_2(i_deriv))
                                    END DO
                                 END DO
                              END DO
                           END DO
                        END IF

                        DEALLOCATE (p_deriv)
                     END DO !lk
                  END DO !li
               END DO !lj

            END DO !kpgf
         END DO !jpgf
      END DO !ipgf

      IF (PRESENT(der_abc_1_ext)) der_abc_1_ext = der_abc_1_ext_prv
      IF (PRESENT(der_abc_2_ext)) der_abc_2_ext = der_abc_2_ext_prv

      DEALLOCATE (params)

   END SUBROUTINE eri_3center_derivs

! **************************************************************************************************
!> \brief Sets the internals of the cp_libint_t object for derivatives of integrals of type (k|ji)
!> \param lib ..
!> \param ri ...
!> \param rj ...
!> \param rk ...
!> \param zeti ...
!> \param zetj ...
!> \param zetk ...
!> \param li_max ...
!> \param lj_max ...
!> \param lk_max ...
!> \param potential_parameter ...
!> \param params_in ...
!> \param params_out ...
!> \note The use of params_in and params_out comes from the fact that one might have to swap
!>       centers 3 and 4 because of angular momenta and pretty much all the parameters of libint
!>       remain the same upon such a change => might avoid recomputing things over and over again
! **************************************************************************************************
   SUBROUTINE set_params_3c_deriv(lib, ri, rj, rk, zeti, zetj, zetk, li_max, lj_max, lk_max, &
                                  potential_parameter, params_in, params_out)

      TYPE(cp_libint_t), INTENT(INOUT)                   :: lib
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: ri, rj, rk
      REAL(dp), INTENT(IN)                               :: zeti, zetj, zetk
      INTEGER, INTENT(IN), OPTIONAL                      :: li_max, lj_max, lk_max
      TYPE(libint_potential_type), INTENT(IN), OPTIONAL  :: potential_parameter
      TYPE(params_3c), OPTIONAL, POINTER                 :: params_in, params_out

      INTEGER                                            :: l
      LOGICAL                                            :: use_gamma
      REAL(dp)                                           :: gammaq, omega2, omega_corr, omega_corr2, &
                                                            prefac, R, S1234, T, tmp
      REAL(dp), ALLOCATABLE, DIMENSION(:)                :: Fm
      TYPE(params_3c), POINTER                           :: params

      IF (PRESENT(params_in)) THEN
         params => params_in

      ELSE
         params => params_out

         params%m_max = li_max + lj_max + lk_max + 1
         gammaq = zeti + zetj
         params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
         params%ZetapEtaInv = 1._dp/(zetk + gammaq)

         params%Q = (zeti*ri + zetj*rj)*params%EtaInv
         params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
         params%Rho = zetk*gammaq/(zetk + gammaq)

         params%Fm = 0.0_dp
         SELECT CASE (potential_parameter%potential_type)
         CASE (do_potential_coulomb)
            T = params%Rho*SUM((params%Q - rk)**2)
            S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
            prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234

            CALL fgamma(params%m_max, T, params%Fm)
            params%Fm = prefac*params%Fm
         CASE (do_potential_truncated)
            R = potential_parameter%cutoff_radius*SQRT(params%Rho)
            T = params%Rho*SUM((params%Q - rk)**2)
            S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
            prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234

            CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
            CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
            IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
            params%Fm = prefac*params%Fm
         CASE (do_potential_short)
            T = params%Rho*SUM((params%Q - rk)**2)
            S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
            prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234

            CALL fgamma(params%m_max, T, params%Fm)

            omega2 = potential_parameter%omega**2
            omega_corr2 = omega2/(omega2 + params%Rho)
            omega_corr = SQRT(omega_corr2)
            T = T*omega_corr2
            ALLOCATE (Fm(prim_data_f_size))

            CALL fgamma(params%m_max, T, Fm)
            tmp = -omega_corr
            DO l = 1, params%m_max + 1
               params%Fm(l) = params%Fm(l) + Fm(l)*tmp
               tmp = tmp*omega_corr2
            END DO
            params%Fm = prefac*params%Fm
         CASE (do_potential_mix_cl_trunc)
            R = potential_parameter%cutoff_radius*SQRT(params%Rho)
            T = params%Rho*SUM((params%Q - rk)**2)
            S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
            prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234

            CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
            CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
            IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)

            ALLOCATE (Fm(prim_data_f_size))
            CALL fgamma(params%m_max, T, Fm)
            DO l = 1, params%m_max + 1
               params%Fm(l) = params%Fm(l) &
                              *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
                              - Fm(l)*potential_parameter%scale_longrange
            END DO
            DEALLOCATE (Fm)

            omega2 = potential_parameter%omega**2
            omega_corr2 = omega2/(omega2 + params%Rho)
            omega_corr = SQRT(omega_corr2)
            T = T*omega_corr2

            ALLOCATE (Fm(prim_data_f_size))
            CALL fgamma(params%m_max, T, Fm)
            tmp = omega_corr
            DO l = 1, params%m_max + 1
               params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
               tmp = tmp*omega_corr2
            END DO
            params%Fm = prefac*params%Fm
         CASE (do_potential_id)
            S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
                        - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
            prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234

            params%Fm(:) = prefac
         CASE DEFAULT
            CPABORT("Requested operator NYI")
         END SELECT

      END IF

      CALL cp_libint_set_params_eri_deriv(lib, rk, rk, rj, ri, rk, &
                                          params%Q, params%W, zetk, 0.0_dp, zetj, zeti, params%ZetaInv, &
                                          params%EtaInv, params%ZetapEtaInv, params%Rho, params%m_max, params%Fm)

   END SUBROUTINE set_params_3c_deriv

! **************************************************************************************************
!> \brief Computes the 2-center electron repulsion integrals (a|b) for a given set of cartesian
!>        gaussian orbitals
!> \param int_ab the integrals as array of cartesian orbitals (allocated before hand)
!> \param la_min ...
!> \param la_max ...
!> \param npgfa ...
!> \param zeta ...
!> \param rpgfa ...
!> \param ra ...
!> \param lb_min ...
!> \param lb_max ...
!> \param npgfb ...
!> \param zetb ...
!> \param rpgfb ...
!> \param rb ...
!> \param dab ...
!> \param lib the libint_t object for evaluation (assume that it is initialized outside)
!> \param potential_parameter the info about the potential
!> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
!>       the libint library must be static initialized, and in case of truncated Coulomb operator,
!>       the latter must be initialized too
! **************************************************************************************************
   SUBROUTINE eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
                          lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
                          dab, lib, potential_parameter)

      REAL(dp), DIMENSION(:, :), INTENT(INOUT)           :: int_ab
      INTEGER, INTENT(IN)                                :: la_min, la_max, npgfa
      REAL(dp), DIMENSION(:), INTENT(IN)                 :: zeta, rpgfa
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: ra
      INTEGER, INTENT(IN)                                :: lb_min, lb_max, npgfb
      REAL(dp), DIMENSION(:), INTENT(IN)                 :: zetb, rpgfb
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: rb
      REAL(dp), INTENT(IN)                               :: dab
      TYPE(cp_libint_t), INTENT(INOUT)                   :: lib
      TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter

      INTEGER                                            :: a_mysize(1), a_offset, a_start, &
                                                            b_offset, b_start, i, ipgf, j, jpgf, &
                                                            li, lj, ncoa, ncob, p1, p2
      REAL(dp)                                           :: dr_ab, zeti, zetj
      REAL(dp), DIMENSION(:), POINTER                    :: p_work

      NULLIFY (p_work)

      dr_ab = 0.0_dp

      IF (potential_parameter%potential_type == do_potential_truncated .OR. &
          potential_parameter%potential_type == do_potential_short .OR. &
          potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
         dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
      ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
         dr_ab = 1000000.0_dp
      END IF

      !Looping over the pgfs
      DO ipgf = 1, npgfa
         zeti = zeta(ipgf)
         a_start = (ipgf - 1)*ncoset(la_max)

         DO jpgf = 1, npgfb
            zetj = zetb(jpgf)
            b_start = (jpgf - 1)*ncoset(lb_max)

            !screening
            IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE

            CALL set_params_2c(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)

            DO li = la_min, la_max
               a_offset = a_start + ncoset(li - 1)
               ncoa = nco(li)
               DO lj = lb_min, lb_max
                  b_offset = b_start + ncoset(lj - 1)
                  ncob = nco(lj)

                  a_mysize(1) = ncoa*ncob
                  CALL cp_libint_get_2eris(li, lj, lib, p_work, a_mysize)

                  DO j = 1, ncob
                     p1 = (j - 1)*ncoa
                     DO i = 1, ncoa
                        p2 = p1 + i
                        int_ab(a_offset + i, b_offset + j) = p_work(p2)
                     END DO
                  END DO

               END DO
            END DO

         END DO
      END DO

   END SUBROUTINE eri_2center

! **************************************************************************************************
!> \brief Sets the internals of the cp_libint_t object for integrals of type (k|j)
!> \param lib ..
!> \param rj ...
!> \param rk ...
!> \param zetj ...
!> \param zetk ...
!> \param lj_max ...
!> \param lk_max ...
!> \param potential_parameter ...
! **************************************************************************************************
   SUBROUTINE set_params_2c(lib, rj, rk, zetj, zetk, lj_max, lk_max, potential_parameter)

      TYPE(cp_libint_t), INTENT(INOUT)                   :: lib
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: rj, rk
      REAL(dp), INTENT(IN)                               :: zetj, zetk
      INTEGER, INTENT(IN)                                :: lj_max, lk_max
      TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter

      INTEGER                                            :: l, op
      LOGICAL                                            :: use_gamma
      REAL(dp)                                           :: omega2, omega_corr, omega_corr2, prefac, &
                                                            R, T, tmp
      REAL(dp), ALLOCATABLE, DIMENSION(:)                :: Fm
      TYPE(params_2c)                                    :: params

      !The internal structure of libint2 is based on 4-center integrals
      !For 2-center, two of those are dummy centers
      !The integral is assumed to be (k|j) where the centers are ordered as:
      !k -> 1, j -> 3 and (the centers #2 & #4 are dummy centers)

      !Note: some variable of 4-center integrals simplify due to dummy centers:
      !      P -> rk, gammap -> zetk
      !      Q -> rj, gammaq -> zetj

      op = potential_parameter%potential_type
      params%m_max = lj_max + lk_max
      params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
      params%ZetapEtaInv = 1._dp/(zetk + zetj)

      params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
      params%Rho = zetk*zetj/(zetk + zetj)

      params%Fm = 0.0_dp
      SELECT CASE (op)
      CASE (do_potential_coulomb)
         T = params%Rho*SUM((rj - rk)**2)
         prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
         CALL fgamma(params%m_max, T, params%Fm)
         params%Fm = prefac*params%Fm
      CASE (do_potential_truncated)
         R = potential_parameter%cutoff_radius*SQRT(params%Rho)
         T = params%Rho*SUM((rj - rk)**2)
         prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)

         CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
         CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
         IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
         params%Fm = prefac*params%Fm
      CASE (do_potential_short)
         T = params%Rho*SUM((rj - rk)**2)
         prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)

         CALL fgamma(params%m_max, T, params%Fm)

         omega2 = potential_parameter%omega**2
         omega_corr2 = omega2/(omega2 + params%Rho)
         omega_corr = SQRT(omega_corr2)
         T = T*omega_corr2
         ALLOCATE (Fm(prim_data_f_size))

         CALL fgamma(params%m_max, T, Fm)
         tmp = -omega_corr
         DO l = 1, params%m_max + 1
            params%Fm(l) = params%Fm(l) + Fm(l)*tmp
            tmp = tmp*omega_corr2
         END DO
         params%Fm = prefac*params%Fm
      CASE (do_potential_mix_cl_trunc)
         R = potential_parameter%cutoff_radius*SQRT(params%Rho)
         T = params%Rho*SUM((rj - rk)**2)
         prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)

         CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
         CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
         IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)

         ALLOCATE (Fm(prim_data_f_size))
         CALL fgamma(params%m_max, T, Fm)
         DO l = 1, params%m_max + 1
            params%Fm(l) = params%Fm(l) &
                           *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
                           - Fm(l)*potential_parameter%scale_longrange
         END DO
         DEALLOCATE (Fm)

         omega2 = potential_parameter%omega**2
         omega_corr2 = omega2/(omega2 + params%Rho)
         omega_corr = SQRT(omega_corr2)
         T = T*omega_corr2

         ALLOCATE (Fm(prim_data_f_size))
         CALL fgamma(params%m_max, T, Fm)
         tmp = omega_corr
         DO l = 1, params%m_max + 1
            params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
            tmp = tmp*omega_corr2
         END DO
         params%Fm = prefac*params%Fm
      CASE (do_potential_id)

         prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
         params%Fm(:) = prefac
      CASE DEFAULT
         CPABORT("Requested operator NYI")
      END SELECT

      CALL cp_libint_set_params_eri(lib, rk, rk, rj, rj, params%ZetaInv, params%EtaInv, &
                                    params%ZetapEtaInv, params%Rho, rk, rj, params%W, &
                                    params%m_max, params%Fm)

   END SUBROUTINE set_params_2c

! **************************************************************************************************
!> \brief Helper function to compare libint_potential_types
!> \param potential1 first potential
!> \param potential2 second potential
!> \return Boolean whether both potentials are equal
! **************************************************************************************************
   PURE FUNCTION compare_potential_types(potential1, potential2) RESULT(equals)
      TYPE(libint_potential_type), INTENT(IN)            :: potential1, potential2
      LOGICAL                                            :: equals

      IF (potential1%potential_type /= potential2%potential_type) THEN
         equals = .FALSE.
      ELSE
         equals = .TRUE.
         SELECT CASE (potential1%potential_type)
         CASE (do_potential_short, do_potential_long)
            IF (potential1%omega /= potential2%omega) equals = .FALSE.
         CASE (do_potential_truncated)
            IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .FALSE.
         CASE (do_potential_mix_cl_trunc)
            IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .FALSE.
            IF (potential1%omega /= potential2%omega) equals = .FALSE.
            IF (potential1%scale_coulomb /= potential2%scale_coulomb) equals = .FALSE.
            IF (potential1%scale_longrange /= potential2%scale_longrange) equals = .FALSE.
         END SELECT
      END IF

   END FUNCTION compare_potential_types

!> \brief Computes the 2-center derivatives of the electron repulsion integrals (a|b) for a given
!>        set of cartesian gaussian orbitals. Returns the derivatives wrt to the first center
!> \param der_ab the derivatives as array of cartesian orbitals (allocated before hand)
!> \param la_min ...
!> \param la_max ...
!> \param npgfa ...
!> \param zeta ...
!> \param rpgfa ...
!> \param ra ...
!> \param lb_min ...
!> \param lb_max ...
!> \param npgfb ...
!> \param zetb ...
!> \param rpgfb ...
!> \param rb ...
!> \param dab ...
!> \param lib the libint_t object for evaluation (assume that it is initialized outside)
!> \param potential_parameter the info about the potential
!> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
!>       the libint library must be static initialized, and in case of truncated Coulomb operator,
!>       the latter must be initialized too
! **************************************************************************************************
   SUBROUTINE eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
                                 lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
                                 dab, lib, potential_parameter)

      REAL(dp), DIMENSION(:, :, :), INTENT(INOUT)        :: der_ab
      INTEGER, INTENT(IN)                                :: la_min, la_max, npgfa
      REAL(dp), DIMENSION(:), INTENT(IN)                 :: zeta, rpgfa
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: ra
      INTEGER, INTENT(IN)                                :: lb_min, lb_max, npgfb
      REAL(dp), DIMENSION(:), INTENT(IN)                 :: zetb, rpgfb
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: rb
      REAL(dp), INTENT(IN)                               :: dab
      TYPE(cp_libint_t), INTENT(INOUT)                   :: lib
      TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter

      INTEGER                                            :: a_mysize(1), a_offset, a_start, &
                                                            b_offset, b_start, i, i_deriv, ipgf, &
                                                            j, jpgf, li, lj, ncoa, ncob, p1, p2
      INTEGER, DIMENSION(3)                              :: permute
      REAL(dp)                                           :: dr_ab, zeti, zetj
      REAL(dp), DIMENSION(:, :), POINTER                 :: p_deriv

      NULLIFY (p_deriv)

      permute = [4, 5, 6]

      dr_ab = 0.0_dp

      IF (potential_parameter%potential_type == do_potential_truncated .OR. &
          potential_parameter%potential_type == do_potential_short .OR. &
          potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
         dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
      ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
         dr_ab = 1000000.0_dp
      END IF

      !Looping over the pgfs
      DO ipgf = 1, npgfa
         zeti = zeta(ipgf)
         a_start = (ipgf - 1)*ncoset(la_max)

         DO jpgf = 1, npgfb
            zetj = zetb(jpgf)
            b_start = (jpgf - 1)*ncoset(lb_max)

            !screening
            IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE

            CALL set_params_2c_deriv(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)

            DO li = la_min, la_max
               a_offset = a_start + ncoset(li - 1)
               ncoa = nco(li)
               DO lj = lb_min, lb_max
                  b_offset = b_start + ncoset(lj - 1)
                  ncob = nco(lj)

                  a_mysize(1) = ncoa*ncob
                  CALL cp_libint_get_2eri_derivs(li, lj, lib, p_deriv, a_mysize)

                  DO i_deriv = 1, 3
                     DO j = 1, ncob
                        p1 = (j - 1)*ncoa
                        DO i = 1, ncoa
                           p2 = p1 + i
                           der_ab(a_offset + i, b_offset + j, i_deriv) = p_deriv(p2, permute(i_deriv))
                        END DO
                     END DO
                  END DO

                  DEALLOCATE (p_deriv)
               END DO
            END DO

         END DO
      END DO

   END SUBROUTINE eri_2center_derivs

! **************************************************************************************************
!> \brief Sets the internals of the cp_libint_t object for derivatives of integrals of type (k|j)
!> \param lib ..
!> \param rj ...
!> \param rk ...
!> \param zetj ...
!> \param zetk ...
!> \param lj_max ...
!> \param lk_max ...
!> \param potential_parameter ...
! **************************************************************************************************
   SUBROUTINE set_params_2c_deriv(lib, rj, rk, zetj, zetk, lj_max, lk_max, potential_parameter)

      TYPE(cp_libint_t), INTENT(INOUT)                   :: lib
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: rj, rk
      REAL(dp), INTENT(IN)                               :: zetj, zetk
      INTEGER, INTENT(IN)                                :: lj_max, lk_max
      TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter

      INTEGER                                            :: l, op
      LOGICAL                                            :: use_gamma
      REAL(dp)                                           :: omega2, omega_corr, omega_corr2, prefac, &
                                                            R, T, tmp
      REAL(dp), ALLOCATABLE, DIMENSION(:)                :: Fm
      TYPE(params_2c)                                    :: params

      !The internal structure of libint2 is based on 4-center integrals
      !For 2-center, two of those are dummy centers
      !The integral is assumed to be (k|j) where the centers are ordered as:
      !k -> 1, j -> 3 and (the centers #2 & #4 are dummy centers)

      !Note: some variable of 4-center integrals simplify due to dummy centers:
      !      P -> rk, gammap -> zetk
      !      Q -> rj, gammaq -> zetj

      op = potential_parameter%potential_type
      params%m_max = lj_max + lk_max + 1
      params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
      params%ZetapEtaInv = 1._dp/(zetk + zetj)

      params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
      params%Rho = zetk*zetj/(zetk + zetj)

      params%Fm = 0.0_dp
      SELECT CASE (op)
      CASE (do_potential_coulomb)
         T = params%Rho*SUM((rj - rk)**2)
         prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
         CALL fgamma(params%m_max, T, params%Fm)
         params%Fm = prefac*params%Fm
      CASE (do_potential_truncated)
         R = potential_parameter%cutoff_radius*SQRT(params%Rho)
         T = params%Rho*SUM((rj - rk)**2)
         prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)

         CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
         CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
         IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
         params%Fm = prefac*params%Fm
      CASE (do_potential_short)
         T = params%Rho*SUM((rj - rk)**2)
         prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)

         CALL fgamma(params%m_max, T, params%Fm)

         omega2 = potential_parameter%omega**2
         omega_corr2 = omega2/(omega2 + params%Rho)
         omega_corr = SQRT(omega_corr2)
         T = T*omega_corr2
         ALLOCATE (Fm(prim_data_f_size))

         CALL fgamma(params%m_max, T, Fm)
         tmp = -omega_corr
         DO l = 1, params%m_max + 1
            params%Fm(l) = params%Fm(l) + Fm(l)*tmp
            tmp = tmp*omega_corr2
         END DO
         params%Fm = prefac*params%Fm
      CASE (do_potential_mix_cl_trunc)
         R = potential_parameter%cutoff_radius*SQRT(params%Rho)
         T = params%Rho*SUM((rj - rk)**2)
         prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)

         CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
         CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
         IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)

         ALLOCATE (Fm(prim_data_f_size))
         CALL fgamma(params%m_max, T, Fm)
         DO l = 1, params%m_max + 1
            params%Fm(l) = params%Fm(l) &
                           *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
                           - Fm(l)*potential_parameter%scale_longrange
         END DO
         DEALLOCATE (Fm)

         omega2 = potential_parameter%omega**2
         omega_corr2 = omega2/(omega2 + params%Rho)
         omega_corr = SQRT(omega_corr2)
         T = T*omega_corr2

         ALLOCATE (Fm(prim_data_f_size))
         CALL fgamma(params%m_max, T, Fm)
         tmp = omega_corr
         DO l = 1, params%m_max + 1
            params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
            tmp = tmp*omega_corr2
         END DO
         params%Fm = prefac*params%Fm
      CASE (do_potential_id)

         prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
         params%Fm(:) = prefac
      CASE DEFAULT
         CPABORT("Requested operator NYI")
      END SELECT

      CALL cp_libint_set_params_eri_deriv(lib, rk, rk, rj, rj, rk, rj, params%W, zetk, 0.0_dp, &
                                          zetj, 0.0_dp, params%ZetaInv, params%EtaInv, &
                                          params%ZetapEtaInv, params%Rho, &
                                          params%m_max, params%Fm)

   END SUBROUTINE set_params_2c_deriv

END MODULE libint_2c_3c

